home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / gamma_inc.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  14.8 KB  |  514 lines
  1. /* specfunc/gamma_inc.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_erf.h>
  26. #include <gsl/gsl_sf_exp.h>
  27. #include <gsl/gsl_sf_log.h>
  28. #include <gsl/gsl_sf_gamma.h>
  29.  
  30. #include "error.h"
  31.  
  32. /* The dominant part,
  33.  * D(a,x) := x^a e^(-x) / Gamma(a+1)
  34.  */
  35. static
  36. int
  37. gamma_inc_D(const double a, const double x, gsl_sf_result * result)
  38. {
  39.   if(a < 10.0) {
  40.     double lnr;
  41.     gsl_sf_result lg;
  42.     gsl_sf_lngamma_e(a+1.0, &lg);
  43.     lnr = a * log(x) - x - lg.val;
  44.     result->val = exp(lnr);
  45.     result->err = 2.0 * GSL_DBL_EPSILON * (fabs(lnr) + 1.0) * fabs(result->val);
  46.     return GSL_SUCCESS;
  47.   }
  48.   else {
  49.     double mu = (x-a)/a;
  50.     double term1;
  51.     gsl_sf_result gstar;
  52.     gsl_sf_result ln_term;
  53.     gsl_sf_log_1plusx_mx_e(mu, &ln_term);  /* log(1+mu) - mu */
  54.     gsl_sf_gammastar_e(a, &gstar);
  55.     term1 = exp(a*ln_term.val)/sqrt(2.0*M_PI*a);
  56.     result->val  = term1/gstar.val;
  57.     result->err  = 2.0 * GSL_DBL_EPSILON * (fabs(a*ln_term.val) + 1.0) * fabs(result->val);
  58.     result->err += gstar.err/fabs(gstar.val) * fabs(result->val);
  59.     return GSL_SUCCESS;
  60.   }
  61. }
  62.  
  63.  
  64. /* P series representation.
  65.  */
  66. static
  67. int
  68. gamma_inc_P_series(const double a, const double x, gsl_sf_result * result)
  69. {
  70.   const int nmax = 5000;
  71.  
  72.   gsl_sf_result D;
  73.   int stat_D = gamma_inc_D(a, x, &D);
  74.  
  75.   double sum  = 1.0;
  76.   double term = 1.0;
  77.   int n;
  78.   for(n=1; n<nmax; n++) {
  79.     term *= x/(a+n);
  80.     sum  += term;
  81.     if(fabs(term/sum) < GSL_DBL_EPSILON) break;
  82.   }
  83.  
  84.   result->val  = D.val * sum;
  85.   result->err  = D.err * fabs(sum);
  86.   result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  87.  
  88.   if(n == nmax)
  89.     GSL_ERROR ("error", GSL_EMAXITER);
  90.   else
  91.     return stat_D;
  92. }
  93.  
  94.  
  95. /* Q large x asymptotic
  96.  */
  97. static
  98. int
  99. gamma_inc_Q_large_x(const double a, const double x, gsl_sf_result * result)
  100. {
  101.   const int nmax = 5000;
  102.  
  103.   gsl_sf_result D;
  104.   const int stat_D = gamma_inc_D(a, x, &D);
  105.  
  106.   double sum  = 1.0;
  107.   double term = 1.0;
  108.   double last = 1.0;
  109.   int n;
  110.   for(n=1; n<nmax; n++) {
  111.     term *= (a-n)/x;
  112.     if(fabs(term/last) > 1.0) break;
  113.     if(fabs(term/sum)  < GSL_DBL_EPSILON) break;
  114.     sum  += term;
  115.     last  = term;
  116.   }
  117.  
  118.   result->val  = D.val * (a/x) * sum;
  119.   result->err  = D.err * fabs((a/x) * sum);
  120.   result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  121.  
  122.   if(n == nmax)
  123.     GSL_ERROR ("error", GSL_EMAXITER);
  124.   else
  125.     return stat_D;
  126. }
  127.  
  128.  
  129. /* Uniform asymptotic for x near a, a and x large.
  130.  * See [Temme, p. 285]
  131.  * FIXME: need c1 coefficient
  132.  */
  133. static
  134. int
  135. gamma_inc_Q_asymp_unif(const double a, const double x, gsl_sf_result * result)
  136. {
  137.   const double rta = sqrt(a);
  138.   const double eps = (x-a)/a;
  139.  
  140.   gsl_sf_result ln_term;
  141.   const int stat_ln = gsl_sf_log_1plusx_mx_e(eps, &ln_term);  /* log(1+eps) - eps */
  142.   const double eta  = eps * sqrt(-2.0*ln_term.val/(eps*eps));
  143.  
  144.   gsl_sf_result erfc;
  145.  
  146.   double R;
  147.   double c0, c1;
  148.  
  149.   gsl_sf_erfc_e(eta*M_SQRT2*rta, &erfc);
  150.  
  151.   if(fabs(eps) < GSL_ROOT5_DBL_EPSILON) {
  152.     c0 = -1.0/3.0 + eps*(1.0/12.0 - eps*(23.0/540.0 - eps*(353.0/12960.0 - eps*589.0/30240.0)));
  153.     c1 = 0.0;
  154.   }
  155.   else {
  156.     double rt_term;
  157.     rt_term = sqrt(-2.0 * ln_term.val/(eps*eps));
  158.     c0 = (1.0 - 1.0/rt_term)/eps;
  159.     c1 = 0.0;
  160.   }
  161.  
  162.   R = exp(-0.5*a*eta*eta)/(M_SQRT2*M_SQRTPI*rta) * (c0 + c1/a);
  163.  
  164.   result->val  = 0.5 * erfc.val + R;
  165.   result->err  = GSL_DBL_EPSILON * fabs(R * 0.5 * a*eta*eta) + 0.5 * erfc.err;
  166.   result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  167.  
  168.   return stat_ln;
  169. }
  170.  
  171.  
  172. /* Continued fraction for Q.
  173.  *
  174.  * Q(a,x) = D(a,x) a/x F(a,x)
  175.  *             1   (1-a)/x  1/x  (2-a)/x   2/x  (3-a)/x
  176.  *  F(a,x) =  ---- ------- ----- -------- ----- -------- ...
  177.  *            1 +   1 +     1 +   1 +      1 +   1 +
  178.  *
  179.  * Hans E. Plesser, 2002-01-22 (hans dot plesser at itf dot nlh dot no):
  180.  *
  181.  * Since the Gautschi equivalent series method for CF evaluation may lead 
  182.  * to singularities, I have replaced it with the modified Lentz algorithm
  183.  * given in
  184.  *
  185.  * I J Thompson and A R Barnett
  186.  * Coulomb and Bessel Functions of Complex Arguments and Order
  187.  * J Computational Physics 64:490-509 (1986)
  188.  *
  189.  * In consequence, gamma_inc_Q_CF_protected() is now obsolete and has been 
  190.  * removed. 
  191.  *
  192.  * Identification of terms between the above equation for F(a, x) and
  193.  * the first equation in the appendix of Thompson&Barnett is as follows:
  194.  *
  195.  *    b_0 = 0, b_n = 1 for all n > 0
  196.  *
  197.  *    a_1 = 1
  198.  *    a_n = (n/2-a)/x    for n even
  199.  *    a_n = (n-1)/(2x)   for n odd
  200.  *
  201.  */
  202.  
  203. static
  204. int
  205. gamma_inc_Q_CF(const double a, const double x, gsl_sf_result * result)
  206. {
  207.   const int    nmax  =  5000;
  208.   const double small =  gsl_pow_3 (GSL_DBL_EPSILON);
  209.  
  210.   gsl_sf_result D;
  211.   const int stat_D = gamma_inc_D(a, x, &D);
  212.  
  213.   double hn = 1.0;           /* convergent */
  214.   double Cn = 1.0 / small;
  215.   double Dn = 1.0;
  216.   int n;
  217.  
  218.   /* n == 1 has a_1, b_1, b_0 independent of a,x, 
  219.      so that has been done by hand                */
  220.   for ( n = 2 ; n < nmax ; n++ ) {
  221.     double an;
  222.     double delta;
  223.  
  224.     if(GSL_IS_ODD(n))
  225.       an = 0.5*(n-1)/x;
  226.     else
  227.       an = (0.5*n-a)/x;
  228.  
  229.     Dn = 1.0 + an * Dn;
  230.     if ( fabs(Dn) < small )
  231.       Dn = small;
  232.     Cn = 1.0 + an/Cn;
  233.     if ( fabs(Cn) < small )
  234.       Cn = small;
  235.     Dn = 1.0 / Dn;
  236.     delta = Cn * Dn;
  237.     hn *= delta;
  238.     if(fabs(delta-1) < GSL_DBL_EPSILON) break;
  239.   }
  240.  
  241.   result->val  = D.val * (a/x) * hn;
  242.   result->err  = D.err * fabs((a/x) * hn);
  243.   result->err += GSL_DBL_EPSILON * (2.0 + 0.5*n) * fabs(result->val);
  244.  
  245.   if(n == nmax)
  246.     GSL_ERROR ("error", GSL_EMAXITER);
  247.   else
  248.     return stat_D;
  249. }
  250.  
  251. /* Useful for small a and x. Handles the subtraction analytically.
  252.  */
  253. static
  254. int
  255. gamma_inc_Q_series(const double a, const double x, gsl_sf_result * result)
  256. {
  257.   double term1;  /* 1 - x^a/Gamma(a+1) */
  258.   double sum;    /* 1 + (a+1)/(a+2)(-x)/2! + (a+1)/(a+3)(-x)^2/3! + ... */
  259.   int stat_sum;
  260.   double term2;  /* a temporary variable used at the end */
  261.  
  262.   {
  263.     /* Evaluate series for 1 - x^a/Gamma(a+1), small a
  264.      */
  265.     const double pg21 = -2.404113806319188570799476;  /* PolyGamma[2,1] */
  266.     const double lnx  = log(x);
  267.     const double el   = M_EULER+lnx;
  268.     const double c1 = -el;
  269.     const double c2 = M_PI*M_PI/12.0 - 0.5*el*el;
  270.     const double c3 = el*(M_PI*M_PI/12.0 - el*el/6.0) + pg21/6.0;
  271.     const double c4 = -0.04166666666666666667
  272.                        * (-1.758243446661483480 + lnx)
  273.                        * (-0.764428657272716373 + lnx)
  274.                * ( 0.723980571623507657 + lnx)
  275.                * ( 4.107554191916823640 + lnx);
  276.     const double c5 = -0.0083333333333333333
  277.                        * (-2.06563396085715900 + lnx)
  278.                * (-1.28459889470864700 + lnx)
  279.                * (-0.27583535756454143 + lnx)
  280.                * ( 1.33677371336239618 + lnx)
  281.                * ( 5.17537282427561550 + lnx);
  282.     const double c6 = -0.0013888888888888889
  283.                        * (-2.30814336454783200 + lnx)
  284.                        * (-1.65846557706987300 + lnx)
  285.                        * (-0.88768082560020400 + lnx)
  286.                        * ( 0.17043847751371778 + lnx)
  287.                        * ( 1.92135970115863890 + lnx)
  288.                        * ( 6.22578557795474900 + lnx);
  289.     const double c7 = -0.00019841269841269841
  290.                        * (-2.5078657901291800 + lnx)
  291.                        * (-1.9478900888958200 + lnx)
  292.                        * (-1.3194837322612730 + lnx)
  293.                        * (-0.5281322700249279 + lnx)
  294.                        * ( 0.5913834939078759 + lnx)
  295.                        * ( 2.4876819633378140 + lnx)
  296.                        * ( 7.2648160783762400 + lnx);
  297.     const double c8 = -0.00002480158730158730
  298.                        * (-2.677341544966400 + lnx)
  299.                        * (-2.182810448271700 + lnx)
  300.                        * (-1.649350342277400 + lnx)
  301.                        * (-1.014099048290790 + lnx)
  302.                        * (-0.191366955370652 + lnx)
  303.                        * ( 0.995403817918724 + lnx)
  304.                        * ( 3.041323283529310 + lnx)
  305.                        * ( 8.295966556941250 + lnx);
  306.     const double c9 = -2.75573192239859e-6
  307.                        * (-2.8243487670469080 + lnx)
  308.                        * (-2.3798494322701120 + lnx)
  309.                        * (-1.9143674728689960 + lnx)
  310.                        * (-1.3814529102920370 + lnx)
  311.                        * (-0.7294312810261694 + lnx)
  312.                        * ( 0.1299079285269565 + lnx)
  313.                        * ( 1.3873333251885240 + lnx)
  314.                        * ( 3.5857258865210760 + lnx)
  315.                        * ( 9.3214237073814600 + lnx);
  316.     const double c10 = -2.75573192239859e-7
  317.                        * (-2.9540329644556910 + lnx)
  318.                        * (-2.5491366926991850 + lnx)
  319.                        * (-2.1348279229279880 + lnx)
  320.                        * (-1.6741881076349450 + lnx)
  321.                        * (-1.1325949616098420 + lnx)
  322.                        * (-0.4590034650618494 + lnx)
  323.                        * ( 0.4399352987435699 + lnx)
  324.                        * ( 1.7702236517651670 + lnx)
  325.                        * ( 4.1231539047474080 + lnx)
  326.                        * ( 10.342627908148680 + lnx);
  327.  
  328.     term1 = a*(c1+a*(c2+a*(c3+a*(c4+a*(c5+a*(c6+a*(c7+a*(c8+a*(c9+a*c10)))))))));
  329.   }
  330.  
  331.   {
  332.     /* Evaluate the sum.
  333.      */
  334.     const int nmax = 5000;
  335.     double t = 1.0;
  336.     int n;
  337.     sum = 1.0;
  338.  
  339.     for(n=1; n<nmax; n++) {
  340.       t *= -x/(n+1.0);
  341.       sum += (a+1.0)/(a+n+1.0)*t;
  342.       if(fabs(t/sum) < GSL_DBL_EPSILON) break;
  343.     }
  344.     
  345.     if(n == nmax)
  346.       stat_sum = GSL_EMAXITER;
  347.     else
  348.       stat_sum = GSL_SUCCESS;
  349.   }
  350.  
  351.   term2 = (1.0 - term1) * a/(a+1.0) * x * sum;
  352.   result->val  = term1 + term2;
  353.   result->err  = GSL_DBL_EPSILON * (fabs(term1) + 2.0*fabs(term2));
  354.   result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  355.   return stat_sum;
  356. }
  357.  
  358.  
  359. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  360.  
  361. int
  362. gsl_sf_gamma_inc_Q_e(const double a, const double x, gsl_sf_result * result)
  363. {
  364.   if(a <= 0.0 || x < 0.0) {
  365.     DOMAIN_ERROR(result);
  366.   }
  367.   else if(x == 0.0) {
  368.     result->val = 1.0;
  369.     result->err = 0.0;
  370.     return GSL_SUCCESS;
  371.   }
  372.   else if(x <= 0.5*a) {
  373.     /* If the series is quick, do that. It is
  374.      * robust and simple.
  375.      */
  376.     gsl_sf_result P;
  377.     int stat_P = gamma_inc_P_series(a, x, &P);
  378.     result->val  = 1.0 - P.val;
  379.     result->err  = P.err;
  380.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  381.     return stat_P;
  382.   }
  383.   else if(a >= 1.0e+06 && (x-a)*(x-a) < a) {
  384.     /* Then try the difficult asymptotic regime.
  385.      * This is the only way to do this region.
  386.      */
  387.     return gamma_inc_Q_asymp_unif(a, x, result);
  388.   }
  389.   else if(a < 0.2 && x < 5.0) {
  390.     /* Cancellations at small a must be handled
  391.      * analytically; x should not be too big
  392.      * either since the series terms grow
  393.      * with x and log(x).
  394.      */
  395.     return gamma_inc_Q_series(a, x, result);
  396.   }
  397.   else if(a <= x) {
  398.     if(x <= 1.0e+06) {
  399.       /* Continued fraction is excellent for x >~ a.
  400.        * We do not let x be too large when x > a since
  401.        * it is somewhat pointless to try this there;
  402.        * the function is rapidly decreasing for
  403.        * x large and x > a, and it will just
  404.        * underflow in that region anyway. We 
  405.        * catch that case in the standard
  406.        * large-x method.
  407.        */
  408.       return gamma_inc_Q_CF(a, x, result);
  409.     }
  410.     else {
  411.       return gamma_inc_Q_large_x(a, x, result);
  412.     }
  413.   }
  414.   else {
  415.     if(a < 0.8*x) {
  416.       /* Continued fraction again. The convergence
  417.        * is a little slower here, but that is fine.
  418.        * We have to trade that off against the slow
  419.        * convergence of the series, which is the
  420.        * only other option.
  421.        */
  422.       return gamma_inc_Q_CF(a, x, result);
  423.     }
  424.     else {
  425.       gsl_sf_result P;
  426.       int stat_P = gamma_inc_P_series(a, x, &P);
  427.       result->val  = 1.0 - P.val;
  428.       result->err  = P.err;
  429.       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  430.       return stat_P;
  431.     }
  432.   }
  433. }
  434.  
  435.  
  436. int
  437. gsl_sf_gamma_inc_P_e(const double a, const double x, gsl_sf_result * result)
  438. {
  439.   if(a <= 0.0 || x < 0.0) {
  440.     DOMAIN_ERROR(result);
  441.   }
  442.   else if(x == 0.0) {
  443.     result->val = 0.0;
  444.     result->err = 0.0;
  445.     return GSL_SUCCESS;
  446.   }
  447.   else if(x < 20.0 || x < 0.5*a) {
  448.     /* Do the easy series cases. Robust and quick.
  449.      */
  450.     return gamma_inc_P_series(a, x, result);
  451.   }
  452.   else if(a > 1.0e+06 && (x-a)*(x-a) < a) {
  453.     /* Crossover region. Note that Q and P are
  454.      * roughly the same order of magnitude here,
  455.      * so the subtraction is stable.
  456.      */
  457.     gsl_sf_result Q;
  458.     int stat_Q = gamma_inc_Q_asymp_unif(a, x, &Q);
  459.     result->val  = 1.0 - Q.val;
  460.     result->err  = Q.err;
  461.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  462.     return stat_Q;
  463.   }
  464.   else if(a <= x) {
  465.     /* Q <~ P in this area, so the
  466.      * subtractions are stable.
  467.      */
  468.     gsl_sf_result Q;
  469.     int stat_Q;
  470.     if(a > 0.2*x) {
  471.       stat_Q = gamma_inc_Q_CF(a, x, &Q);
  472.     }
  473.     else {
  474.       stat_Q = gamma_inc_Q_large_x(a, x, &Q);
  475.     }
  476.     result->val  = 1.0 - Q.val;
  477.     result->err  = Q.err;
  478.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  479.     return stat_Q;
  480.   }
  481.   else {
  482.     if((x-a)*(x-a) < a) {
  483.       /* This condition is meant to insure
  484.        * that Q is not very close to 1,
  485.        * so the subtraction is stable.
  486.        */
  487.       gsl_sf_result Q;
  488.       int stat_Q = gamma_inc_Q_CF(a, x, &Q);
  489.       result->val  = 1.0 - Q.val;
  490.       result->err  = Q.err;
  491.       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  492.       return stat_Q;
  493.     }
  494.     else {
  495.       return gamma_inc_P_series(a, x, result);
  496.     }
  497.   }
  498. }
  499.  
  500.  
  501. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  502.  
  503. #include "eval.h"
  504.  
  505. double gsl_sf_gamma_inc_P(const double a, const double x)
  506. {
  507.   EVAL_RESULT(gsl_sf_gamma_inc_P_e(a, x, &result));
  508. }
  509.  
  510. double gsl_sf_gamma_inc_Q(const double a, const double x)
  511. {
  512.   EVAL_RESULT(gsl_sf_gamma_inc_Q_e(a, x, &result));
  513. }
  514.